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Abstract: Since the characterization of the eye's monochromatic aberration 
fluctuations in 2001, the power spectrum has remained the most widely used 
method for analyzing their dynamics. However, the power spectrum does 
not capture the complexities of the fluctuations. We measured the 
monochromatic aberration dynamics of six subjects using a Shack- 
Hartmann sensor sampling at 21 Hz. We characterized the dynamics using 
techniques from chaos theory. We found that the attractor embedding 
dimension for all aberrations, for all subjects, was equal to three. The 
embedding lag averaged across aberrations and subjects was 0.31 ± 0.07 s. 
The Lyapunov exponent of the rms wavefront error was positive for each 
subject, with an average value of 0.44 ± 0.15 (j.m/s. This indicates that the 
aberration dynamics are chaotic. Implications for future modeling are 
discussed. 
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1. Introduction 

The dynamics of the monochromatic aberrations of the human eye, during steady-state 
fixation, were first characterized just over a decade ago [1]. In this study, power spectrum 
analysis revealed that the fluctuations can be described by an inverse power law. Since this 
time, power spectrum based measures, in particular the slope of the power spectrum, have 
remained the most widely used methods for characterizing monochromatic aberration 
dynamics; see for example [2-9]. However, characterizing the dynamics in this way is likely 
to be an oversimplification. For example, it has recently been demonstrated that the dynamics 
are multifractal and characterized by a range of statistically distributed exponents [10]. This 
has implications for modeling aberration dynamics, and determining differences both between 
subjects, and for the same subject with differing experimental conditions. 

In the investigations of other physiological signals, such as the heartbeat, more advanced 
analytical methods are commonplace. One such technique is chaos theory [11]. Chaos theory 
characterizes nonlinear dynamical systems that are sensitive to initial conditions [12]. Unlike 
the usual meaning of the word chaos, chaotic systems have underlying laws that determine 
their behavior. The reason they appear to be 'chaotic' is because of their sensitivity to initial 
conditions — the so called Butterfly Effect. Figure 1 shows an example of a chaotic time 
series. Examples of chaotic signals in physiology are the heartbeat [13], brain waves [14], 
pupillary unrest [15], and blood pressure variations [16]. 



1 




Fig. 1. Example of a chaotic time series. The generating equation is the so-called logistic 
equation: x t + ; = kx,(l-x,), which describes the relative size of a population over time, t. In the 
graphs the growth rate, k, is 3.8. A change in the initial relative population, x h by only 0.1% 
results in a divergence of the plots. 

The aim of this work was to determine if chaos is present in the fluctuations of the eye's 
monochromatic aberrations. 

2. Method 

2.1. Subjects 

The aberrations of the right eye of six subjects were measured in this study. Aberration 
measurements were performed over the natural pupil size of each subject. All subjects gave 
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informed consent to take part in the study. Their demographics are shown in Table 1 . Those 
who required spectacles wore them during the experiment. 



Table 1. Subject demographics 



Subject 


Age (yrs) 


Pupil (mm) 


Refraction (RE) 


KH 


29 


5.5 


-1.75/-0.50 x 90 


EM 


34 


5.0 


-6.00/-0.50 x 90 


JC 


32 


6.0 


Piano 


YP 


24 


5.0 


Piano 


CS 


2? 


4.5 


-6.00/-0.50 x 60 


cv 


28 


4.0 


-4.25/- 1.25 x 5 



2.2. Aberration measurements 

The aberrations were measured over a 23 s time period using an open-view Shack-Hartmann 
sensor sampling at 21 Hz. The aberrations of subjects KH and CV were measured over 11.5 s 
as they found it difficult to maintain fixation on the stimulus for 23 s. The subjects were 
stabilized using a bite bar. A simplified schematic of the sensor is shown in Fig. 2. For further 
details of the system see [17]. The stimulus was a Maltese cross displayed on an LCD monitor 
at an accommodative demand of 0.4 D. It subtended 11.32 minutes of arc at the eye and had a 
luminance of 255cd/m 2 . The left eye of each subject remained open and fixated on the 
stimulus. For each subject the time evolution of each aberration (up to and including 5th radial 
order excluding tip and tilt), and the rms wavefront error were calculated. Blinks were 
removed using a cubic spline function [5]. The analysis was carried out using custom written 
software in Matlab. 




Fig. 2. Shack-Hartmann sensor for measurement of the eye's aberrations. 

3. Chaos theory analysis 

A feature of chaos is sensitivity to initial conditions. This manifests itself as the exponential 
divergence (separation) of neighboring trajectories in phase space [12]. Phase space is a plot 
in which each axis is effectively a variable needed to specify the 'state' of the system over time 
[12]. The plotted points map out the so-called attractor. Figure 3(a) shows a schematic of the 
phase space plot for the logistic equation in Fig. 1. Figure 3(b) shows a schematic of the 
change in separation of neighboring trajectories versus time for such a chaotic time series. For 
a chaotic series there is an initial increase in the rate of divergence for a limited period of 
time. This is the limit of predictability [12]. Following this, the divergence no longer increases 
as the trajectories have diverged so much they are effectively moving randomly with respect 
to each other [12]. The slope of the linear region is the Lyapunov exponent. For a stochastic 
time series the divergence plot is flat and the Lyapunov exponent is zero, and for a non- 
chaotic time series in which trajectories converge, the exponent is negative. The aim of the 
analysis was to reconstruct the attractor of the aberration fluctuations for each subject and 
determine the Lyapunov exponents. 
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Fig. 3. (a) Schematic of the phase space plot for the logistic equation in Fig. 1. The divergence 
of neighboring trajectories is also shown, (b) A schematic of the change in the separation of 
neighboring trajectories versus time for such a chaotic time series. 

3. 1. Phase space reconstruction 

The first step in chaos theory analysis is reconstruction of the phase space. To reconstruct the 
phase space from the one dimensional time series, time -delay embedding was used [18]. In 
this method, a point in J-dimensional phase space is given by 



p(0 = [x(t), x(t + r),...,x(t + (d- l)r)], 



(1) 



where the embedding lag r = t 0 + iAt, with At being the time between frames [18]. The 
embedding dimension d is equal to the number of axes and number of 'sub-series'. The 
number of data points, n, per sub-series is given by 



N-(d-V)xi, 



(2) 



where N is the total number of data points, and i is the embedding lag in units of data points. 
3.1.1. Embedding lag 

The lag, r, must be large enough so that the data points in each sub-series are not correlated, 
otherwise there is redundancy, as one sub-series can predict another. In this case the data 
points will accumulate along a diagonal line [18]. However, the lag must not be so large that 
the data points in each sub-series are completely unrelated [18]. In this case, the phase space 
plot will look like random scatter [12]. The optimum lag was determined by using the first 
minimum of the mutual information [11]. For two time series x and y, the mutual information 
is given by 



(3) 



where Pxy is the joint probability, and P x and P Y are the marginal probabilities [12]. For the 
purpose of determining r, the two time series are the original series and the lagged series. 
Hence 
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\- / u, ( P(x(t),x(t + r)) 
I(t)= > P(x(0,x(? + z-))log 2 — v — v — 



x(f),*(f+0 



(4) 



A high value indicates that the two series contain the same information, i.e. the data points 
in one series can predict that of the other. The mutual information method is preferred over 
using an autocorrelation-based method as autocorrelation relies on a linear basis [12]. 

For each subject, for each Zernike aberration, I(t) was calculated for a lag of i = 1 to 30 
data points. Prior to calculating I(z), the time course signals were also detrended as trends 
cause unwanted correlations [12]. 

3.1.2. Embedding dimension 

The embedding dimension, d, was calculated using a false nearest neighbors (FNN) approach 
[19]. This involves measuring how the distances between neighboring points in phase space 
vary with increasing dimension. Suppose a data point located in d-dimensional space is given 

by 

p(0 =[x(t),x(t + T),...,x(t + (d-l)r)], (5) 
and it's nearest neighbor is given by 

PNN ( { NN )=l X ( t NNl X ( t NN+ T l-> X ( { NN + (A ~ l ) T )] ■ ( 6 ) 

The distance between the points is 

K (p,p w ) = Z w x(t NN + (k - \)t )f . (7) 

If the embedding dimension is too small the points will appear to be closer than they 
actually are, as illustrated in Fig. 4(a). Hence as the embedding dimension increases the 
measured distance will increase until the correct embedding dimension is found (Fig. 4(b)), 
where after the distance remains constant (Fig. 4(c)). The separation of points when the 
dimension increases by one is given by 

r m (P> Pnn ) = R ] (P> Pnn ) + [ x (t + dt) - x(t m +dr)f. (8) 
And so the change in distance can be calculated as 



< d(P.PNw)-^(P.P«w) _ \ x (t + dr)-x(t m +dT)\ 

V ^ 2 (p,Pnn) ^(p.Pnn) 

A point is considered a FNN if two conditions are satisfied: 



(9) 



^ 2 + i(P>Pnn)--^(P>Pnn) >r ^ ( 10 ) 



and 

where R A is the size of the attractor, which can be estimated using the standard deviation of 
the time series [19]. This second condition accounts for the fact that the nearest neighbor of a 
point may not necessarily be close to it. Without this extra condition, the number of false 
nearest neighbors starts to rise again as d increases beyond the appropriate dimension. 
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For a data series, the nearest neighbor of each point is determined for a given dimension. 
For each pair of points, Eqs. (10) and (11) are then evaluated. The embedding dimension is 
then determined from a plot of the percentage of false nearest neighbors versus dimension. 
For each subject, for each Zernike aberration, the percentage of false nearest neighbors was 
determined for d = 1:10. A value of 15 was used for R Tot [20]. A Tot was set to 5. This was the 
lowest value such that the percentage of false nearest neighbors did not start to rise again as d 
continued to increase. 




Fig. 4. Principle of the false nearest neighbors method to determine the embedding dimension. 
The correct embedding dimension in the illustration is 2. A point and its neighbor are separated 
by a distance Rr n ,e- (a) Data points are incorrectly embedded in one-dimensional phase space. 
The points appear to be closer to each other than they actually are {RMeas < Rrrue), and therefore 
the nearest neighbor is false, (b) Data points are correctly embedded in two-dimensional phase 
space and the true distance is determined, (c) Data points embedded in a higher dimension. 
Again the measured distance is equal to the true distance. 

3.2. Lyapunov exponent 

The number of Lyapunov exponents is equal to the number of axes in the phase space 
reconstruction [12]. However in practice, only the largest exponent is calculated as this 
dominates the trajectory divergence [21]. There are several algorithms which calculate the 
largest Lyapunov exponent from the reconstructed attractor [12]. We used the algorithm of 
Rosenstein et al. as it is robust to noise and suitable for small data sets [21]. The divergence of 
trajectories at a time t is given by 

d(t) = Ce L ', (12) 

where C is a constant and L is the largest Lyapunov exponent. For each data point in the 
reconstructed attractor, its nearest neighbor is found. The nearest neighbor has to be separated 
in time by more than the average period of the time series. This constraint allows one to 
consider the two points to represent two different nearby trajectories. The mean period is 
calculated as the inverse of the mean frequency of the power spectrum [21]. From Eq. (12), 
for the j th pair of neighbors, and time step i, the separation is 

d j (i) = Ce L(m , (13) 

where At is the time between frames. Hence 

\nd j {i) = \nC + L(iM). (14) 

For each aberration and the rms wavefront error for each subject, the left-hand side of Eq. 
(14) was evaluated for i = 1:60. The Lyapunov exponent was determined from the slope of the 
linear rise, if applicable, in the divergence averaged across neighbors, versus time (iAt) [21]. 
To determine the region of the linear rise, the data were fitted using two straight lines given 

by 

y x = m l x(l:n) + C 1 (15) 

and 
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y 2 = m 2 x(n : end) + C 2 , 



(16) 



where n is the break-point time step, n was determined as the value which minimizes 

m 2 x(n) + C 2 -(m l x(n) + C l ). (17) 
Hence m 2 is the Lyapunov exponent, and (nAf) is the limit of predictability. 
4. Results 

Figure 5 shows the mutual information for each subject averaged across Zernike coefficients. 
The mean lag across all subjects was 0.31 ± 0.07 s. The average percentage of false nearest 
neighbors averaged across all aberrations for each subject is shown in Fig. 6. The embedding 
dimension was taken as the lag at which the number of false nearest neighbors was < 5%. This 
gave an embedding dimension of 3 for each subject. The embedding dimension for each 
individual aberration for each subject was also 3. Figure 7 shows the phase space 
reconstruction for each subject for the rms wavefront error. 



KH EM JC 




0 10 20 30 0 10 20 30 0 10 20 30 



Lag Lag Lag 

Fig. 5. The mutual information averaged across Zernike coefficients for each subject. The lag is 
in units of data points. * Indicates the location of the first minimum, and hence the embedding 
lag. 
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Fig. 6. Graph of the percentage of FNN averaged across all Zernike coefficients for each 
subject. * Indicates the embedding dimension, which was taken as the point where the lag was 
<5%. 

KH EM JC 
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Figure 8 shows the average divergence across the attractor for the rms wavefront error for 
each subject. Only the first 40 data points, corresponding to a time span of up to 2 s, have 
been included for clarity. The linear rise is also shown. For comparison purposes the 
divergence of the rms wavefront error, randomly reordered in time, is also plotted. This so- 
called random-shuffle surrogate data preserves the statistical properties of the data and is used 
to confirm that the data is not random but chaotic [14]. As can be seen from the plots, 
reordering the data eliminates the linear rise, and the divergence plot is characteristic of 
stochastic data. The mean limit of predictability for the rms wavefront error was 0.99 ± 0.20 s, 
and the mean Lyapunov exponent was 0.44 ± 0.15 um/s. The lag, limit of predictability, and 
Lyapunov exponent for each subject for the rms wavefront error are shown in Table 2. 



KH EM JC 




0 0.5 1 1.5 2 0 0.5 1 1.5 2 0 0.5 1 1.5 2 

Time (s) Time (s) Time (s) 



YP CS CV 




0 0.5 1 1.5 2 0 0.5 1 1.5 2 0 0.5 1 1.5 2 

Time (s) Time (s) Time (s) 



Fig. 8. Average divergence of neighboring trajectories for the rms wavefront error for each 
subject. The linear-rise region, limit of predictability (dotted line), and divergence of random- 
shuffle surrogate data, are also shown. The random-shuffle data has been displaced by +1 um/s 
for each subject for clarity. 



Table 2. Lag, limit of predictability P, m „ and Lyapunov exponent L rms , for the rms 
wavefront error of each subject, and averaged across subjects 



Subject 


Lag (s) 


Prms (S) 


Lms (um/s) 


KH 


0.24 


0.82 


0.51 


EM 


0.34 


1.02 


0.27 


JC 


0.34 


1.16 


0.35 


YP 


0.39 


1.21 


0.52 


CS 


0.34 


1.02 


0.33 


CV 


0.19 


0.68 


0.66 


Mean 


0.31 ±0.07 


0.99 ± 0.20 


0.44 ±0.15 



Figure 9 shows the limit of predictability and Lyapunov exponent for each individual 
Zernike aberration coefficient averaged across subjects. The lags used in the phase space 
reconstruction for each subject were obtained from Fig. 5, the first minimum of the average 
mutual information, and the embedding dimension was 3. The attractors for individual 
Zernike aberrations for each subject were similar to those shown in Fig. 7. 
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Fig. 9. (a) The limit of predictability for each individual Zernike aberration coefficient 
averaged across subjects, (b) Corresponding Lyapunov exponents. The lags used in the phase 
space reconstruction for each subject were obtained from Fig. 5, the minimum of the average 
mutual information. Error bars show ± S.D. 

5. Discussion 

In this study we measured the monochromatic aberration dynamics of 6 subjects during 
steady-state fixation. We used techniques from chaos theory to reconstruct the attractors and 
determine the Lyapunov exponents. 

5. 1. Chaotic nature of aberration dynamics and implications for modeling 

There are currently very few models of ocular aberration dynamics [22]. Modeling of 
aberration dynamics is important, for example, in developing algorithms to better control 
closed-loop adaptive optics systems [23]. From Fig. 8 it can be seen that in comparison to the 
randomly shuffled data, all subjects show an initial rise in the divergence of neighboring 
trajectories for the rms wavefront error. The mean Lyapunov exponent was 0.44 ± 0.15 |im/s. 
All individual Zernike aberration coefficients also had a positive Lyapunov exponent. This 
indicates that like other physiological signals, the fluctuations in the aberrations are chaotic 
[11,13-16]. Consequently, the aberration fluctuations are not stochastic and there are 
underlying equations that determine their behavior. The phase space portraits from the rms 
wavefront error shown in Fig. 7 are similar to those of other physiological systems and are 
characteristic of a so-called strange attractor [12]. 

The mean lag averaged across aberrations and subjects for the rms wavefront error was 
0.31 ± 0.07 s. This value indicates the memory of the system generating the fluctuations, as 
beyond this time the data are reasonably uncorrelated. The average limit of predictability for 
the rms wavefront error was 0.99 ± 0.20 s. Hence after this time period it is difficult to 
quantify the separation of trajectories in phase space. The limit of predictability was similar 
across all individual aberration coefficients. For all subjects, the embedding dimension for all 
Zernike aberration coefficients was 3. This indicates that there are ~3 variables controlling the 
fluctuations [20]. Hence, when developing models, at least three variables should be included. 

5.2 Possible inputs to the chaotic dynamics 

5.2.1. The heartbeat 

We speculate that one of the major inputs into the chaotic dynamics of the eyes' aberrations is 
the effect of the heartbeat. A number of studies have shown that the aberration dynamics show 
some correlation with the heartbeat [4,8,24], and there is much evidence that the heartbeat is 
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itself a chaotic signal, see for example [13,20]. The heartbeat has a wide-spread impact on the 
eye, and its effect on aberrations is likely to mediate itself via a number of mechanisms. These 
include fundus pulsation [25], corneal pulsation [26], and the effect on the crystalline lens 
owing to blood flow through the ciliary body [24]. Furthermore, the heartbeat is correlated 
with fixational eye movements [27,28], changes in intraocular pressure [29,30], and 
fluctuations in pupil size [31]. These events are likely to alter the shape of the eye and affect 
the aberrations. Although these events show some correlation with the heartbeat, we note that 
their dynamic properties unrelated to the heartbeat may also affect the chaos observed in this 
study. Chaos has been observed in fixational eye movements [32], blood pressure control 
[16], and pupillary hippus [15]. 

The average values of the Lyapunov exponents found here are similar to that of the 
heartbeat [20]. Future work will involve measuring simultaneously the heartbeat signal 
alongside the ocular aberrations, and processing both signals using chaos theory. The likely 
inpact of the heartbeat on the chaotic nature of the aberration dynamics raises the interesting 
possibility that changes in the chaos in the heartbeat, for example due to ill health [33,34], 
may reveal itself as changes in the chaotic properties of aberration dynamics. 

5.2.2. Accommodative microfluctuations 

Another plausible major input into the chaotic dynamics is the microfluctuations in 
accommodation [35]. Aside from the impact of the heartbeat [36], unlike other aberrations, 
there is evidence that components of accommodation fluctuations are actively changed by the 
accommodation control system, see for example [37]. Specifically, the magnitude of the 
fluctuations increases in response to a decrease in blur sensitivity, and so are affected by 
retinal image quality [38-41]. Several authors have proposed that the accommodation control 
system monitors the rate of change of image contrast with microfluctuations in defocus, to 
help the eye stay in focus on a stationary target, and also to determine the response to changes 
in object distance [42]. 

In this study we also analyzed the microfluctuations in accommodation using the defocus 
and spherical aberration Zernike coefficients [43]. We found the attractor to be similar to that 
of the rms wavefront error. The mean limit of predictability across subjects was 1.16 ± 0.35 s, 
and the mean Lyapunov exponent was 0.36 ± 0.26 D/s. As the other aberrations in the eye 
show some correlation with accommodation microfluctuations, this may partly explain the 
chaos observed in their fluctuations [44,45]. From Fig. 9(b) the Lyapunov exponent of 
Zernike defocus is notably larger than that of the other aberrations. Owing to the variability 
between subjects, this difference was not statistically significant however. The larger value 
may be due to the active input of the accommodation control system to this aberration. Future 
work will involve altering retinal image quality, with adaptive optics for example [38], and 
determining how the chaotic properties of microfluctuations in accommodation and other 
aberrations change. 

5.2.3. Tear film fluctuations 

Changes in tear film are known to have some impact on the aberration dynamics, see for 
example [46]. To examine the effect of tear film fluctuations we fitted subject EM with a 
scleral contact lens [47]. The front of the lens was dried once the lens was in place on the eye. 
This lens preserves the tear film between the eye and the lens. Figure 10(a) shows the 
fluctuations in the rms wavefront error with and without the lens in place. For this subject, 
there were no blinks during the measurement without and with the scleral lens. Hence we 
were comparing the effect of the tear film breaking up. Figure 10(b) shows the corresponding 
divergence plots. The scleral lens did not affect the limit of predictability. The Lyapunov 
exponent was 0.33 |im/s for the scleral lens case, in comparison to a value of 0.27 (xm/s for 
the no lens case. The increase of the Lyapunov exponent for the scleral lens case indicates that 
the fluctuations in tear film may be a source of noise that 'dilutes' the chaos, as a flat 
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Fig. 10. (a) Time course of the rms wavefront error for subject EM with and without a scleral 
lens. Both signals have been detrended and the plot for the no lens case has been displaced by 
an amount of +0.15 urn for clarity, (b) Corresponding divergence plots. Again the no lens case 
has been displaced for clarity (+0.55 um/s). 

divergence plot indicates noise. Hence the Lyapunov exponents obtained in this study should 
be regarded as a lower limit. 

5.2.4. Others 

Another possible input into the results is cognitive demand. It has been shown that cognitive 
demand affects both the heartbeat and accommodation system [48]. Future work will include 
analyzing the chaotic properties of aberration dynamics with different cognitive tasks. 

5.3. Potential impact of instrument noise, eye movement artifacts and blink removal artifacts 

The Shack-Hartmann measurements are not noise free. There are noise sources intrinsic to the 
instrument itself, such as read noise, and sources external to the instrument owing to 
movement of the eye relative to the system. The effect of noise is to shorten the length of time 
of the linear rising region in the divergence plots [21], and also eventually flatten the 
divergence plots. Hence, like the Lyapunov exponents, the limits of predictability obtained in 
this study should be considered as a lower bound of the noise free (true) limits of 
predictability. 

To estimate the intrinsic noise level, measurements were taken on an artificial eye 
consisting of a 20 mm focal length lens and a flat matte surface acting as a retina [17]. 
Measurements were performed over a 5 mm pupil diameter and the laser power was lowered 
such that the amount of light returning from the artificial eye was the same as that returning 
from the human eye. Figure 11(a) shows the time course of the rms wavefront error for both 
the artificial eye and subject YP (whose measurements were also recorded over a 5 mm 
pupil). As can be seen from the plot, the fluctuations in the artificial eye are considerably 
smaller than that of YP. Hence the contribution of intrinsic noise to the observed results is 
likely to be negligible. The signal to noise, measured as the ratio of the standard deviation of 
the rms wavefront error for YP to that of the artificial eye is 13. Figure 11(b) shows the 
divergence plots for both signals, which confirms that the results for the artificial eye are not 
chaotic but merely stochastic. 

Fixational eye movements cause an effective translation of the pupil relative to the lenslet 
array [9]. Sahin has found that the mean pupil translation of healthy subjects during a 50 ms 
time span is around 40 ± 10 (xm [49]. To determine the relative effect of eye movements, for 
each Shack-Hartmann measurement for subject KH, we calculated the rms wavefront error 
with the pupil shifted by +40 (xm and -40 |im. Subject KH was chosen as they had the 
smallest fluctuations in their rms wavefront error and so pupil movements will have had a 
larger relative effect. Figure 12(a) shows the rms wavefront error time course for the original 
time series, assuming the eye's pupil was static, and the resulting plots for a pupil movement 
of +40 (xm and -40 |im. As can be seen from the plot the differences owing to pupil 
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Fig. 11. (a) Time course of the rms wavefront error for subject YP and an artificial eye. Both 
signals have been shifted to 0 urn for comparison purposes, (b) Corresponding divergence 
plots. The artificial eye data has been shifted by +3.1 um/s for clarity. 
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Fig. 12. (a) Time course of the rms wavefront error for subject KH assuming a static pupil and 
the pupil shifted by plus and minus 40 um. (b) Corresponding divergence plots. Plots have 
been displaced for clarity. The +40 um and -40 um results have been shifted by a value of -0.8 
um/s and —1.5 um/s respectively. 

movement are extremely small in comparison to the measured fluctuations assuming a fixed 
pupil. Figure 12(b) shows the divergence plots. The limit of predictability ranged from 0.78 s 
to 0.97 s, and the Lyapunov exponent ranged from 0.50 (im/s to 0.51 |im/s. In comparison to 
the fixed pupil results (P rms = 0.82 s and L rms = 0.51 Lim/s), we therefore estimate the impact of 
pupil movements to be 23% for the limit of predictability and 2% for the Lyapunov exponent. 
Future work will involve simultaneous measurement of the pupil position during the Shack- 
Hartmann recordings. 

Prior to analyzing the data using chaos theory, blinks were removed and the data points 
were replaced using cubic spline interpolation. Typically there were two to three blinks during 
a 23 s time period and one to two blinks during a 11.5 s time period. This amounts to up to 
around only 4% of data points being replaced by a cubic spline function. Hence it is unlikely 
that removal of blinks had a significant impact on the results. 

5.4 Other applications of chaos theory analysis to aberration dynamics 

A change in the chaotic properties of other physiological signals can indicate changes in the 
state of the system, for example owing to pathology [33,34]. Hence chaos theory may have 
application in detecting subtle changes in the aberration dynamics both between subjects, and 
within subjects for different experimental conditions. In particular, chaos theory may have 
utility in the study of micro fluctuations of accommodation and myopia [41]. For example, in 
late-onset myopes, the fluctuations in accommodation display large low frequency 
fluctuations [50]. Methods used to control chaos in the brain and heart could also be applied to 
the eye to restore normal functioning of the accommodation system in such subjects [51,52]. 
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6. Conclusion 

We used techniques from chaos theory to analyze the monochromatic aberration dynamics of 
six subjects during steady-state viewing. We found that the fluctuations in the aberrations are 
chaotic and that the dynamics can be reconstructed in phase space in 3 dimensions; this has 
implications for future modeling. A further application of chaos theory is in studying and 
manipulating the accommodation system, particularly in myopic subjects. 
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